GLMM Part2

Preparations

Load the necessary libraries

library(car) # for regression diagnostics
library(broom) # for tidy output
library(ggfortify) # for model diagnostics
library(sjPlot) # for outputs
library(knitr) # for kable
library(effects) # for partial effects plots
library(ggeffects) # for partial effects plots
library(emmeans) # for estimating marginal means
library(MASS) # for glm.nb
library(MuMIn) # for AICc
library(tidyverse) # for data wrangling
library(nlme)
library(lme4)
library(glmmTMB)
library(broom.mixed)
library(glmmTMB) # for glmmTMB
library(DHARMa) # for residuals and diagnostics
library(performance) # for diagnostic plots
library(see) # for diagnostic plots
library(patchwork) # for multiplot layouts

Scenario

To investigate differential metabolic plasticity in barramundi (Lates calcarifer), Norin, Malte, and Clark (2015) exposed juvenile barramundi to various environmental changes (increased temperature, decreased salinity and increased hypoxia) as well as control conditions. Metabolic plasticity was calculated as the percentage difference in standard metabolic rate between the various treatment conditions and the standard metabolic rate under control conditions. They were interested in whether there was a relationship between metabolic plasticity and typical (control) metabolism and how the different treatment conditions impact on this relationship.

A total of 60 barramundi juveniles were subject to each of the three conditions (high temperature, low salinity and hypoxia) in addition to control conditions. Fish mass was also recorded as a covariate as this is known to influence metabolic parameters.

Barramundi

Sampling design

Format of norin.csv data files

FISHID	MASS	TRIAL	SMR_contr	CHANGE
1	35.69	LowSalinity	5.85	-31.92
2	33.84	LowSalinity	6.53	2.52
3	37.78	LowSalinity	5.66	-6.28
..	..	..	..	..
1	36.80	HighTemperature	5.85	18.32
2	34.98	HighTemperature	6.53	19.06
3	38.38	HighTemperature	5.66	19.03
..	..	..	..	..
1	45.06	Hypoxia	5.85	-18.61
2	43.51	Hypoxia	6.53	-5.37
3	45.11	Hypoxia	5.66	-13.95

FISHID	Categorical listing of the individual fish that are repeatedly sampled
MASS	Mass (g) of barramundi. Covariate in analysis
TRIAL	Categorical listing of the trial (LowSalinity: 10ppt salinity; HighTemperature: 35 degrees; Hypoxia: 45% air-sat. oxygen.
SMR_contr	Standard metabolic rate (mg/h/39.4 g of fish) under control trial conditions (35 ppt salinity, 29 degrees, normoxia)
CHANGE	Percentage difference in Standard metabolic rate (mg/h/39.4 g of fish) between Trial conditions and control adjusted for 'regression to the mean'.

Read in the data

norin <- read_csv("../public/data/norin.csv", trim_ws = TRUE)
glimpse(norin)

## Rows: 180
## Columns: 5
## $ FISHID    <dbl> 1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13, 14, 15, 16, 17, 1…
## $ MASS      <dbl> 35.69, 33.84, 37.78, 26.58, 37.62, 37.68, 30.62, 50.37, 24.9…
## $ TRIAL     <chr> "LowSalinity", "LowSalinity", "LowSalinity", "LowSalinity", …
## $ SMR_contr <dbl> 5.847466, 6.530707, 5.659556, 6.278200, 4.407336, 4.818589, …
## $ CHANGE    <dbl> -31.919389, 2.520929, -6.284968, -4.346675, -3.071329, -15.0…

Since the response variable is a difference between the SMR under control conditions compared to the specific trial conditions (either Low salinity, High temperature or Hypoxia), the only possible distribution we can use is a Gaussian. Well to clarify this, we could use a t-distribution, but this is not really available for frequentist mixed effects routines in R.

Exploratory data analysis

Model formula: \[ y_i \sim{} \mathcal{N}(\mu_i, \sigma^2)\\ \mu_i =\boldsymbol{\beta} \bf{X_i} + \boldsymbol{\gamma} \bf{Z_i} \]

where \(\boldsymbol{\beta}\) and \(\boldsymbol{\gamma}\) are vectors of the fixed and random effects parameters respectively and \(\bf{X}\) is the model matrix representing the overall intercept and effects of temperature and (centred) mean fish size on SDA peak. \(\bf{Z}\) represents a cell means model matrix for the random intercepts associated with individual fish.

Lets begin by declaring the categorical predictors and random effects as factors.

norin <- norin %>% mutate(
  FISHID = factor(FISHID),
  TRIAL = factor(TRIAL)
)

Each of the fish were exposed to each of the three trials (treatments). So lets explore the distributions of responses within each of these trials.

ggplot(norin, aes(y = CHANGE, x = TRIAL)) +
  geom_boxplot()

Conclusions:

• these seem reasonable enough

The researchers considered that physiological plasticity might be effected by the fishes basal metabolic rate. For example, a fish with relatively high metabolism might have less scope to increase this metabolism than a fish with a lower metabolism. Therefore, the SMR under control conditions (just prior to the specific trial) could be added as a continuous covariate.

Doing so would introduce two additional assumptions:

1. linearity: that the relationship between the response and control SMR is linear (if we intend to model it as a linear trend)
2. homogeneity of slopes: that the rate of change in response associated with a one unit change in control SMR is similar for each trial (unless we include an interaction)

ggplot(norin, aes(y = CHANGE, x = SMR_contr, shape = TRIAL, color = TRIAL)) +
  geom_smooth(method = "lm") +
  geom_point()

ggplot(norin, aes(y = CHANGE, x = SMR_contr, shape = TRIAL, color = TRIAL)) +
  geom_smooth() +
  geom_point()

It might also be worth exploring how consistent the trial effect is across the different fish. This can give us an idea of whether the addition of a random intercept/slope model would be useful.

ggplot(norin, aes(y = CHANGE, x = as.numeric(FISHID), color = TRIAL)) +
  geom_point() +
  geom_line()

Conclusions:

• generally, the difference between control and trial SMR is greatest for the High temperature trial.
• there is less consistency about the relative effects of Hypoxia and Low Salinity.
• hence a random intercept/slope model might be useful.

Finally, as an acknowledgement that the metabolic response might be influenced

by the mass of the fish, the researchers contemplated including fish Mass as a continuous covariate. There are numerous alternative ways to incorporate covariates that are known to impact on a response.

• add them as an offset. This will include the covariate in the model, yet will not estimate a partial slope for this predictor. Rather the slope is assumed to be 1. Since the parameter does not need to be estimated, adding a predictor as an offset does not consume any degrees of freedom. However, it does assume that there is a slope of one. This approach can be useful as an alternative to using the response divided by the covariate as the response. That is, it is equivalent to using a response that has been standardised by the covariate, yet still allowing the modelling distribution that would be appropriate for the original response.
• add them as a regular covariate. This will include the covariate in the model and estimate a partial slope parameter just like any other term. Although it will incur a degrees of freedom penalty, it does not assume a slope of 1 and could be modelled as non-linear if appropriate.

Lets explore the relationship between the response and Mass separately for each Trial.

ggplot(norin, aes(y = CHANGE, x = MASS, color = TRIAL)) +
  geom_point() +
  geom_smooth(method = "lm")

Conclusions:

• there does not seem to be much of a relationship between the response and Mass. Nevertheless, in the presence of control SMN, it might still be a useful covariate at explaining some of the unexplained variability (and thus increasing the power of the model).
• there is no evidence that the response and Mass are related in a non-linear manner.

Fit the model

For each of the following modelling alternatives, we will:

• establish whether a random intercept/slope model is useful (based on AICc). This step involves evaluating alternative random effects structures. When doing so, we should use Restricted Maximum Likelihood (REML)
• establish whether MASS is a useful covariate (based on AICc). This step involves evaluating alternative fixed effects structures. When doing so, we must use Maximum Likelihood (ML) rather than Residual Maximum Likelihood (REML).

lme (nlme)

random structures

## Start with random intercept model
norin.lme1 <- lme(CHANGE ~ 1 + offset(MASS),
  random = ~ 1 | FISHID,
  data = norin,
  method = "REML"
)
## Then random intercept and slope model
norin.lme2 <- lme(CHANGE ~ 1 + offset(MASS),
  random = ~ TRIAL | FISHID,
  data = norin,
  method = "REML"
)
AICc(norin.lme1, norin.lme2)

anova(norin.lme1, norin.lme2)

Conclusions:

• the AICc of the random intercept/slope model is lowest
• therefore we will proceed with a random intercept/slope model

fixed structures

## Compare models that estimate partial slope for MASS vs an offset for MASS
## must use ML to compare models that vary in fixed effects
norin.lme2a <- update(norin.lme2,
  . ~ TRIAL * scale(SMR_contr, scale = FALSE) + offset(MASS),
  method = "ML"
)
## The following does not converge
## norin.lme2b <- update(norin.lme2,
##                       .~TRIAL*scale(SMR_contr, scale=FALSE) + MASS,
##                       method = 'ML')
## Use BFGS from optim instead
norin.lme2b <- update(norin.lme2,
  . ~ TRIAL * scale(SMR_contr, scale = FALSE) + MASS,
  method = "ML",
  control = lmeControl(
    opt = "optim",
    optimMethod = "BFGS"
  )
)
norin.lme2c <- update(norin.lme2,
  . ~ TRIAL * scale(SMR_contr, scale = FALSE),
  method = "ML"
)

AICc(norin.lme2a, norin.lme2b, norin.lme2c)

anova(norin.lme2a, norin.lme2b, norin.lme2c)

## Now that we have decided on the structure, we need to ensure that the
## model has been fit with REML
norin.lme2c <- update(norin.lme2c, method = "REML")

Conclusions:

• on the basis of AICc, the model without MASS is considered the most parsimonious (lowest AICc)
• we will proceed with the fixed structure of model 2c

lmer (lme4)

random structure

## Start with random intercept model
norin.lmer1 <- lmer(CHANGE ~ 1 + offset(MASS) + (1 | FISHID),
  data = norin,
  REML = TRUE
)
## Then random intercept and slope model
norin.lmer2 <- lmer(CHANGE ~ 1 + offset(MASS) + (TRIAL | FISHID),
  data = norin,
  REML = TRUE,
  control = lmerControl(check.nobs.vs.nRE = "ignore")
)
AICc(norin.lmer1, norin.lmer2)

anova(norin.lmer1, norin.lmer2)

Conclusions:

• the AICc of the random intercept/slope model is lowest
• therefore we will proceed with a random intercept/slope model

fixed structure

## Compare models that estimate partial slope for MASS vs an offset for MASS
## must use ML to compare models that vary in fixed effects
norin.lmer2a <- update(norin.lmer2,
  . ~ TRIAL * scale(SMR_contr, scale = FALSE) + offset(MASS)
    + (TRIAL | FISHID),
  data = norin,
  REML = FALSE
)
norin.lmer2b <- update(norin.lmer2,
  . ~ TRIAL * scale(SMR_contr, scale = FALSE) + MASS
    + (TRIAL | FISHID),
  data = norin,
  REML = FALSE
)
norin.lmer2c <- update(norin.lmer2,
  . ~ TRIAL * scale(SMR_contr, scale = FALSE)
    + (TRIAL | FISHID),
  data = norin,
  REML = FALSE
)
## Compare these models via AICc
AICc(norin.lmer2a, norin.lmer2b, norin.lmer2c)

## Alternatively,  we can use sequential Likelihood Ratio Tests (LRT)
anova(norin.lmer2a, norin.lmer2b, norin.lmer2c)

## Now that we have decided on the structure, we need to ensure that the
## model has been fit with REML
norin.lmer2c <- update(norin.lmer2c, REML = TRUE)

Conclusions:

• on the basis of AICc, the model without MASS is considered the most parsimonious (lowest AICc)
• we will proceed with the fixed structure of model 2c

glmmTMB

random structure

## Start with random intercept model
norin.glmmTMB1 <- glmmTMB(CHANGE ~ 1 + offset(MASS) + (1 | FISHID),
  data = norin,
  REML = TRUE
)
## Then random intercept and slope model
norin.glmmTMB2 <- glmmTMB(CHANGE ~ 1 + offset(MASS) + (TRIAL | FISHID),
  data = norin,
  REML = TRUE
)
## norin.glmmTMB2 <- glmmTMB(CHANGE ~ 1 + offset(MASS) + (TRIAL|FISHID),
##                           data = as.data.frame(norin),
##                           REML = TRUE,
##                           control=glmmTMBControl(optimizer=optim,
##                                                  optArgs=list(method='BFGS')
##                                                  )
##                           )
norin.glmmTMB2 <- glmmTMB(CHANGE ~ 1 + offset(scale(MASS, scale = FALSE)) + (TRIAL | FISHID),
  data = norin,
  REML = TRUE
)
AICc(norin.glmmTMB1, norin.glmmTMB2)

anova(norin.glmmTMB1, norin.glmmTMB2)

Conclusions:

• the AICc of the random intercept/slope model is lowest
• therefore we will proceed with a random intercept/slope model

fixed structure

## Compare models that estimate partial slope for MASS vs an offset for MASS
norin.glmmTMB2a <- glmmTMB(CHANGE ~ TRIAL * scale(SMR_contr, scale = FALSE) +
  offset(MASS) +
  (TRIAL | FISHID),
data = norin,
REML = FALSE
)
## OR with update
norin.glmmTMB2a <- update(norin.glmmTMB2,
  . ~ TRIAL * scale(SMR_contr, scale = FALSE) +
    offset(scale(MASS, scale = FALSE)) +
    (TRIAL | FISHID),
  REML = FALSE
)
## now with MASS as a partial slope
norin.glmmTMB2b <- glmmTMB(CHANGE ~ TRIAL * scale(SMR_contr, scale = FALSE) +
  scale(MASS, scale = FALSE) +
  (TRIAL | FISHID),
data = norin,
REML = FALSE
)
## OR with update
norin.glmmTMB2b <- update(norin.glmmTMB2,
  . ~ TRIAL * scale(SMR_contr, scale = FALSE) +
    scale(MASS, scale = FALSE) +
    (TRIAL | FISHID),
  REML = FALSE
)
## Finally without MASS
norin.glmmTMB2c <- glmmTMB(CHANGE ~
TRIAL * scale(SMR_contr, scale = FALSE) +
  (TRIAL | FISHID),
data = norin,
REML = FALSE
)
## I have explored a range of alternative optimizers etc and
## none of them resolve the issue.
## The current model will be good enough to explore AICc and it turns out that the model
## converges fine once we run with REML

## Note the following does resolve it (not sure why), yet causes downstream issues
## norin.glmmTMB2c <- glmmTMB(CHANGE~
##                                TRIAL*scale(SMR_contr, scale=FALSE) +
##                               (TRIAL|FISHID),
##                           data = norin,
##                           control=glmmTMBControl(profile=TRUE),
##                          REML=FALSE)
## ## OR with update
## norin.glmmTMB2c <- update(norin.glmmTMB2,
##                           .~TRIAL*scale(SMR_contr, scale=FALSE) +
##                               (TRIAL|FISHID),
##                           control=glmmTMBControl(profile=TRUE),
##                          REML=FALSE)
## Compare these models via AICc
AICc(norin.glmmTMB2a, norin.glmmTMB2b, norin.glmmTMB2c)

## Now that we have decided on the structure, we need to ensure that the
## model has been fit with REML
## It turns out that when the model withough MASS is run using REML, the
## response type residuals are tiny - not sure why.
## We will instead go with the offset(MASS) model (2a)
norin.glmmTMB2c <- update(norin.glmmTMB2a, REML = TRUE)

Conclusions:

• on the basis of AICc, the model without MASS is considered the most parsimonious (lowest AICc)
• we will proceed with the fixed structure of model 2c

Model validation

lme (nlme)

Conclusions:

• whilst the residual plot is OK, the Q-Q plot deviates substantially from a straight line. This suggests that there is more mass in the tails of the distribution than a Gaussian would expect. Although the model is reasonably robust to this, it would be nice to have the option to model these data against a flatter distribution such as a Student t. Unfortunately, this is not available.

plot_model

plot_model(norin.lme2c, type = "diag")[-2] %>% plot_grid()

Conclusions:

• whilst the residual plot is OK, the Q-Q plot deviates substantially from a straight line. This suggests that there is more mass in the tails of the distribution than a Gaussian would expect. Although the model is reasonably robust to this, it would be nice to have the option to model these data against a flatter distribution such as a Student t. Unfortunately, this is not available.

Performance model checking

norin.lme2c %>% performance::check_model()

## Could not compute standard errors from random effects for diagnostic plot.
## Homogeneity of variance could not be computed. Cannot extract residual variance from objects of class 'lme'.

DHARMa residuals

Unfortunately, this did not work…

## norin.resid = simulateResiduals(norin.lme2c,  plot=TRUE)

lmer (lme4)

plot_model

plot_model(norin.lmer2c, type = "diag")[-2] %>% plot_grid()

Conclusions:

• whilst the residual plot is OK, the Q-Q plot deviates substantially from a straight line. This suggests that there is more mass in the tails of the distribution than a Gaussian would expect. Although the model is reasonably robust to this, it would be nice to have the option to model these data against a flatter distribution such as a Student t. Unfortunately, this is not available.

Performance model checking

norin.lmer2c %>% performance::check_model()

Conclusions:

• whilst the residual plot is OK, the Q-Q plot deviates substantially from a straight line. This suggests that there is more mass in the tails of the distribution than a Gaussian would expect. Although the model is reasonably robust to this, it would be nice to have the option to model these data against a flatter distribution such as a Student t. Unfortunately, this is not available.

DHARMa residuals

norin.resid <- norin.lmer2c %>% simulateResiduals(plot = TRUE)

Conclusions:

• in contrast to the regular Q-Q normal plots, the Q-Q normal plots generated from the simulated residuals look fine.
• the residual plot from the simulated residuals look fine
• the diagnostic tests reveal no issues.

glmmTBM (glmmTMB)

plot_model

plot_model(norin.glmmTMB2c, type = "diag")[-2] %>% plot_grid()

Conclusions:

• whilst the residual plot is OK, the Q-Q plot deviates substantially from a straight line. This suggests that there is more mass in the tails of the distribution than a Gaussian would expect. Although the model is reasonably robust to this, it would be nice to have the option to model these data against a flatter distribution such as a Student t. Unfortunately, this is not available.

Performance model checking

norin.glmmTMB2c %>% performance::check_model()

Conclusions:

• whilst the residual plot is OK, the Q-Q plot deviates substantially from a straight line. This suggests that there is more mass in the tails of the distribution than a Gaussian would expect. Although the model is reasonably robust to this, it would be nice to have the option to model these data against a flatter distribution such as a Student t. Unfortunately, this is not available.

DHARMa residuals

norin.resid <- norin.glmmTMB2c %>% simulateResiduals(plot = TRUE)

Conclusions:

• in contrast to the regular Q-Q normal plots, the Q-Q normal plots generated from the simulated residuals look fine.
• the residual plot from the simulated residuals look fine
• the diagnostic tests reveal no issues.

Partial plots model

lme (nlme)

plot_model

norin.lme2c %>% plot_model(type = "eff", terms = c("SMR_contr", "TRIAL"), show.data = TRUE)

We can also explore the fixed effects (as a caterpillar plot).

norin.lme2c %>% plot_model(type = "est")

And similarly, the individual random effects (intercepts of each fish).

# plot_model(norin.lme3a,  type='re')

allEffects

norin.lme2c %>%
  allEffects() %>%
  plot(multiline = TRUE, ci.style = "bands")

## Unfortunately, the use of scale() seems to prevent partial.residuals being calculated
## norin.lme2c %>% allEffects() %>% plot(multiline=TRUE,  ci.style='bands', partial.residuals=TRUE)

ggpredict

norin.lme2c %>%
  ggpredict(c("SMR_contr", "TRIAL")) %>%
  plot(add.data = TRUE)

ggemmeans

norin.lme2c %>%
  ggemmeans(~ SMR_contr * TRIAL) %>%
  plot(add.data = TRUE)

lmer (lmer)

plot_model

norin.lmer2c %>%
  plot_model(type = "eff", show.data = TRUE) %>%
  plot_grid()

norin.lmer2c %>% plot_model(type = "eff", terms = c("SMR_contr", "TRIAL"), show.data = TRUE)

We can also explore the fixed effects (as a caterpillar plot).

norin.lmer2c %>% plot_model(type = "est")

And similarly, the individual random effects (intercepts of each fish).

norin.lmer2c %>% plot_model(type = "re")

allEffects

norin.lmer2c %>%
  allEffects() %>%
  plot(multiline = TRUE, ci.style = "bands")

ggpredict

norin.lmer2c %>%
  ggpredict(c("SMR_contr", "TRIAL")) %>%
  plot(add.data = TRUE)

ggemmeans

norin.lmer2c %>%
  ggemmeans(~ SMR_contr * TRIAL) %>%
  plot(add.data = TRUE)

glmmTMB (glmmTMB)

plot_model

norin.glmmTMB2c %>%
  plot_model(type = "eff", show.data = TRUE) %>%
  plot_grid()

norin.glmmTMB2c %>%
  plot_model(type = "eff", terms = c("SMR_contr", "TRIAL"), show.data = TRUE)

We can also explore the fixed effects (as a caterpillar plot).

norin.glmmTMB2c %>% plot_model(type = "est")

And similarly, the individual random effects (intercepts of each fish).

norin.glmmTMB2c %>% plot_model(type = "re")

allEffects

norin.glmmTMB2c %>%
  allEffects() %>%
  plot(multiline = TRUE, ci.style = "bands")

ggpredict

norin.glmmTMB2c %>%
  ggpredict(c("SMR_contr", "TRIAL")) %>%
  plot(add.data = TRUE)

ggemmeans

norin.glmmTMB2c %>%
  ggemmeans(~ SMR_contr * TRIAL) %>%
  plot(add.data = TRUE)

Model investigation / hypothesis testing

lme (nlme)

summary

norin.lme2c %>% summary()

## Linear mixed-effects model fit by REML
##   Data: norin 
##        AIC      BIC    logLik
##   1593.961 1635.028 -783.9803
## 
## Random effects:
##  Formula: ~TRIAL | FISHID
##  Structure: General positive-definite, Log-Cholesky parametrization
##                  StdDev    Corr         
## (Intercept)      16.323848 (Intr) TRIALH
## TRIALHypoxia     14.926864 -0.601       
## TRIALLowSalinity 30.304028 -0.111 -0.215
## Residual          8.659992              
## 
## Fixed effects:  CHANGE ~ TRIAL + scale(SMR_contr, scale = FALSE) + TRIAL:scale(SMR_contr,      scale = FALSE) 
##                                                      Value Std.Error  DF
## (Intercept)                                       50.38244  2.385594 116
## TRIALHypoxia                                     -50.22621  2.492663 116
## TRIALLowSalinity                                 -42.22318  4.219647 116
## scale(SMR_contr, scale = FALSE)                  -21.55293  3.820913  58
## TRIALHypoxia:scale(SMR_contr, scale = FALSE)      13.51107  3.992402 116
## TRIALLowSalinity:scale(SMR_contr, scale = FALSE)  10.18925  6.758444 116
##                                                     t-value p-value
## (Intercept)                                       21.119453  0.0000
## TRIALHypoxia                                     -20.149618  0.0000
## TRIALLowSalinity                                 -10.006330  0.0000
## scale(SMR_contr, scale = FALSE)                   -5.640779  0.0000
## TRIALHypoxia:scale(SMR_contr, scale = FALSE)       3.384195  0.0010
## TRIALLowSalinity:scale(SMR_contr, scale = FALSE)   1.507633  0.1344
##  Correlation: 
##                                                  (Intr) TRIALH TRIALL s(Ss=F
## TRIALHypoxia                                     -0.620                     
## TRIALLowSalinity                                 -0.215 -0.035              
## scale(SMR_contr, scale = FALSE)                   0.000  0.000  0.000       
## TRIALHypoxia:scale(SMR_contr, scale = FALSE)      0.000  0.000  0.000 -0.620
## TRIALLowSalinity:scale(SMR_contr, scale = FALSE)  0.000  0.000  0.000 -0.215
##                                                  TRIALHs=F
## TRIALHypoxia                                              
## TRIALLowSalinity                                          
## scale(SMR_contr, scale = FALSE)                           
## TRIALHypoxia:scale(SMR_contr, scale = FALSE)              
## TRIALLowSalinity:scale(SMR_contr, scale = FALSE) -0.035   
## 
## Standardized Within-Group Residuals:
##         Min          Q1         Med          Q3         Max 
## -1.66585731 -0.30692682  0.02265975  0.29410064  1.72690964 
## 
## Number of Observations: 180
## Number of Groups: 60

Conclusions:

• the average change in SMR for the High Temperature trial at the average control SMR is 50.38.
• on average, the change in SMR during Hypoxia is 50.23 less at the average control SMR than during the High Temperature trial.
• on average, the change in SMR during Low Salinity is 42.22 less at the average control SMR than during the High Temperature trial.
• during the High Temperature trial, for every one unit change in control SMR, the change in SMR declined by 21.55
• the rate of change in SMR during the Hypoxia trial was 13.51 more than that associated with High Temperature.
• the rate of change in SMR during the Low Salinity trial was 10.19 more than that associated with High Temperature (although this was not found to be significant).
• there is approximately twice as much variation within Fish (replicate level) as there is between Fish (8.66 vs
  16.32).

confint

norin.lme2c %>% intervals()

## Approximate 95% confidence intervals
## 
##  Fixed effects:
##                                                       lower      est.     upper
## (Intercept)                                       45.657469  50.38244  55.10741
## TRIALHypoxia                                     -55.163245 -50.22621 -45.28918
## TRIALLowSalinity                                 -50.580718 -42.22318 -33.86563
## scale(SMR_contr, scale = FALSE)                  -29.201316 -21.55293 -13.90454
## TRIALHypoxia:scale(SMR_contr, scale = FALSE)       5.603613  13.51107  21.41852
## TRIALLowSalinity:scale(SMR_contr, scale = FALSE)  -3.196697  10.18925  23.57520
## 
##  Random Effects:
##   Level: FISHID 
##                                         lower       est.       upper
## sd((Intercept))                     5.0780215 16.3238483  52.4747730
## sd(TRIALHypoxia)                    0.9252388 14.9268645 240.8149035
## sd(TRIALLowSalinity)               15.1035225 30.3040284  60.8026464
## cor((Intercept),TRIALHypoxia)      -0.8169111 -0.6005325  -0.2359666
## cor((Intercept),TRIALLowSalinity)  -0.8197871 -0.1109852   0.7321192
## cor(TRIALHypoxia,TRIALLowSalinity) -0.9832297 -0.2148081   0.9603211
## 
##  Within-group standard error:
##       lower        est.       upper 
##   0.1424309   8.6599921 526.5391714

tidy

norin.lme2c %>% tidy(conf.int = TRUE)

norin.lme2c %>%
  tidy(conf.int = TRUE) %>%
  kable()

effect	group	term	estimate	std.error	df	statistic	p.value	conf.low	conf.high
fixed	fixed	(Intercept)	50.3824384	2.385594	116	21.119453	0.0000000	45.6574691	55.1074076
fixed	fixed	TRIALHypoxia	-50.2262117	2.492663	116	-20.149618	0.0000000	-55.1632454	-45.2891780
fixed	fixed	TRIALLowSalinity	-42.2231759	4.219647	116	-10.006330	0.0000000	-50.5807177	-33.8656340
fixed	fixed	scale(SMR_contr, scale = FALSE)	-21.5529273	3.820913	58	-5.640779	0.0000005	-29.2013157	-13.9045389
fixed	fixed	TRIALHypoxia:scale(SMR_contr, scale = FALSE)	13.5110683	3.992402	116	3.384195	0.0009745	5.6036133	21.4185233
fixed	fixed	TRIALLowSalinity:scale(SMR_contr, scale = FALSE)	10.1892534	6.758444	116	1.507633	0.1343671	-3.1966967	23.5752035
ran_pars	FISHID	sd_(Intercept)	16.3238483	NA	NA	NA	NA	5.0780215	52.4747730
ran_pars	FISHID	cor_TRIALHypoxia.(Intercept)	-0.6005325	NA	NA	NA	NA	-0.8169111	-0.2359666
ran_pars	FISHID	cor_TRIALLowSalinity.(Intercept)	-0.1109852	NA	NA	NA	NA	-0.8197871	0.7321192
ran_pars	FISHID	sd_TRIALHypoxia	14.9268645	NA	NA	NA	NA	0.9252388	240.8149035
ran_pars	FISHID	cor_TRIALLowSalinity.TRIALHypoxia	-0.2148081	NA	NA	NA	NA	-0.9832297	0.9603211
ran_pars	FISHID	sd_TRIALLowSalinity	30.3040284	NA	NA	NA	NA	15.1035225	60.8026464
ran_pars	Residual	sd_Observation	8.6599921	NA	NA	NA	NA	NA	NA

Conclusions:

• the average change in SMR for the High Temperature trial at the average control SMR is 50.38.
• on average, the change in SMR during Hypoxia is 50.23 less at the average control SMR than during the High Temperature trial.
• on average, the change in SMR during Low Salinity is 42.22 less at the average control SMR than during the High Temperature trial.
• during the High Temperature trial, for every one unit change in control SMR, the change in SMR declined by 21.55
• the rate of change in SMR during the Hypoxia trial was 13.51 more than that associated with High Temperature.
• the rate of change in SMR during the Low Salinity trial was 10.19 more than that associated with High Temperature (although this was not found to be significant).
• there is approximately twice as much variation within Fish (replicate level) as there is between Fish (8.66 vs
  16.32).

tab_model

# warning this is only appropriate for html output
norin.lme2c %>% sjPlot::tab_model(show.se = TRUE, show.aic = TRUE)

	CHANGE
Predictors	Estimates	std. Error	CI	p
(Intercept)	50.38	2.39	45.66 – 55.11	<0.001
TRIAL [Hypoxia]	-50.23	2.49	-55.16 – -45.29	<0.001
TRIAL [LowSalinity]	-42.22	4.22	-50.58 – -33.87	<0.001
SMR contr	-21.55	3.82	-29.20 – -13.90	<0.001
TRIAL [Hypoxia] * SMR contr	13.51	3.99	5.60 – 21.42	0.001
TRIAL [LowSalinity] * SMR contr	10.19	6.76	-3.20 – 23.58	0.134
Random Effects
σ2	75.00
τ00 FISHID	266.47
τ11 FISHID.TRIALHypoxia	222.81
τ11 FISHID.TRIALLowSalinity	918.33
ρ01	-0.60
	-0.11
ICC	0.87
N FISHID	60
Observations	180
Marginal R2 / Conditional R2	0.494 / 0.935
AIC	1593.961

lmer (lme4)

summary

norin.lmer2c %>% summary()

## Linear mixed model fit by REML ['lmerMod']
## Formula: CHANGE ~ TRIAL + scale(SMR_contr, scale = FALSE) + (TRIAL | FISHID) +  
##     TRIAL:scale(SMR_contr, scale = FALSE)
##    Data: norin
## Control: lmerControl(check.nobs.vs.nRE = "ignore")
## 
## REML criterion at convergence: 1568
## 
## Scaled residuals: 
##      Min       1Q   Median       3Q      Max 
## -2.30188 -0.42411  0.03131  0.40638  2.38624 
## 
## Random effects:
##  Groups   Name             Variance Std.Dev. Corr       
##  FISHID   (Intercept)      198.27   14.081              
##           TRIALHypoxia      86.41    9.296   -0.60      
##           TRIALLowSalinity 781.94   27.963    0.03 -0.64
##  Residual                  143.19   11.966              
## Number of obs: 180, groups:  FISHID, 60
## 
## Fixed effects:
##                                                  Estimate Std. Error t value
## (Intercept)                                        50.382      2.386  21.119
## TRIALHypoxia                                      -50.226      2.493 -20.150
## TRIALLowSalinity                                  -42.223      4.220 -10.006
## scale(SMR_contr, scale = FALSE)                   -21.553      3.821  -5.641
## TRIALHypoxia:scale(SMR_contr, scale = FALSE)       13.511      3.992   3.384
## TRIALLowSalinity:scale(SMR_contr, scale = FALSE)   10.189      6.758   1.508
## 
## Correlation of Fixed Effects:
##             (Intr) TRIALH TRIALL s(Ss=F TRIALHs=F
## TRIALHypoxi -0.620                               
## TRIALLwSlnt -0.215 -0.035                        
## s(SMR_,s=FA  0.000  0.000  0.000                 
## TRIALH:(s=F  0.000  0.000  0.000 -0.620          
## TRIALLS:s=F  0.000  0.000  0.000 -0.215 -0.035

Conclusions:

• the average change in SMR for the High Temperature trial at the average control SMR is 50.38.
• on average, the change in SMR during Hypoxia is 50.23 less at the average control SMR than during the High Temperature trial.
• on average, the change in SMR during Low Salinity is 42.22 less at the average control SMR than during the High Temperature trial.
• during the High Temperature trial, for every one unit change in control SMR, the change in SMR declined by 21.55
• the rate of change in SMR during the Hypoxia trial was 13.51 more than that associated with High Temperature.
• the rate of change in SMR during the Low Salinity trial was 10.19 more than that associated with High Temperature (although this was not found to be significant).
• there is approximately twice as much variation within Fish (replicate level) as there is between Fish (11.97 vs
  14.08).

confint

norin.lmer2c %>% confint()

##                                                       2.5 %    97.5 %
## .sig01                                             0.000000  21.38850
## .sig02                                                   NA        NA
## .sig03                                                   NA        NA
## .sig04                                                   NA        NA
## .sig05                                                   NA        NA
## .sig06                                                   NA        NA
## .sigma                                                   NA        NA
## (Intercept)                                       45.710775  55.05408
## TRIALHypoxia                                     -55.107544 -45.34488
## TRIALLowSalinity                                 -50.486428 -33.95992
## scale(SMR_contr, scale = FALSE)                  -29.035350 -14.07051
## TRIALHypoxia:scale(SMR_contr, scale = FALSE)       5.692828  21.32930
## TRIALLowSalinity:scale(SMR_contr, scale = FALSE)  -3.045679  23.42418

tidy

norin.lmer2c %>% tidy(conf.int = TRUE)

norin.lmer2c %>%
  tidy(conf.int = TRUE) %>%
  kable()

effect	group	term	estimate	std.error	statistic	conf.low	conf.high
fixed	NA	(Intercept)	50.3824384	2.385590	21.119486	45.706768	55.05811
fixed	NA	TRIALHypoxia	-50.2262117	2.492656	-20.149672	-55.111729	-45.34069
fixed	NA	TRIALLowSalinity	-42.2231759	4.219658	-10.006302	-50.493554	-33.95280
fixed	NA	scale(SMR_contr, scale = FALSE)	-21.5529273	3.820907	-5.640788	-29.041768	-14.06409
fixed	NA	TRIALHypoxia:scale(SMR_contr, scale = FALSE)	13.5110683	3.992391	3.384205	5.686126	21.33601
fixed	NA	TRIALLowSalinity:scale(SMR_contr, scale = FALSE)	10.1892534	6.758463	1.507629	-3.057090	23.43560
ran_pars	FISHID	sd__(Intercept)	14.0807964	NA	NA	NA	NA
ran_pars	FISHID	cor__(Intercept).TRIALHypoxia	-0.5968891	NA	NA	NA	NA
ran_pars	FISHID	cor__(Intercept).TRIALLowSalinity	0.0337670	NA	NA	NA	NA
ran_pars	FISHID	sd__TRIALHypoxia	9.2958585	NA	NA	NA	NA
ran_pars	FISHID	cor__TRIALHypoxia.TRIALLowSalinity	-0.6361548	NA	NA	NA	NA
ran_pars	FISHID	sd__TRIALLowSalinity	27.9632565	NA	NA	NA	NA
ran_pars	Residual	sd__Observation	11.9663529	NA	NA	NA	NA

Conclusions:

• the average change in SMR for the High Temperature trial at the average control SMR is 50.38.
• on average, the change in SMR during Hypoxia is 50.23 less at the average control SMR than during the High Temperature trial.
• on average, the change in SMR during Low Salinity is 42.22 less at the average control SMR than during the High Temperature trial.
• during the High Temperature trial, for every one unit change in control SMR, the change in SMR declined by 21.55
• the rate of change in SMR during the Hypoxia trial was 13.51 more than that associated with High Temperature.
• the rate of change in SMR during the Low Salinity trial was 10.19 more than that associated with High Temperature (although this was not found to be significant).
• there is approximately twice as much variation within Fish (replicate level) as there is between Fish (11.97 vs
  14.08).

tab_model

# warning this is only appropriate for html output
norin.lmer2c %>% sjPlot::tab_model(show.se = TRUE, show.aic = TRUE)

	CHANGE
Predictors	Estimates	std. Error	CI	p
(Intercept)	50.38	2.39	45.67 – 55.09	<0.001
TRIAL [Hypoxia]	-50.23	2.49	-55.15 – -45.31	<0.001
TRIAL [LowSalinity]	-42.22	4.22	-50.55 – -33.89	<0.001
SMR contr	-21.55	3.82	-29.10 – -14.01	<0.001
TRIAL [Hypoxia] * SMR contr	13.51	3.99	5.63 – 21.39	0.001
TRIAL [LowSalinity] * SMR contr	10.19	6.76	-3.15 – 23.53	0.134
Random Effects
σ2	143.19
τ00 FISHID	198.27
τ11 FISHID.TRIALHypoxia	86.41
τ11 FISHID.TRIALLowSalinity	781.94
ρ01	-0.60
	0.03
ICC	0.76
N FISHID	60
Observations	180
Marginal R2 / Conditional R2	0.494 / 0.877
AIC	1593.961

glmmTMB (glmmTMB)

summary

norin.glmmTMB2c %>% summary()

##  Family: gaussian  ( identity )
## Formula:          
## CHANGE ~ TRIAL + scale(SMR_contr, scale = FALSE) + (TRIAL | FISHID) +  
##     TRIAL:scale(SMR_contr, scale = FALSE) + offset(scale(MASS,  
##     scale = FALSE))
## Data: norin
## 
##      AIC      BIC   logLik deviance df.resid 
##   1606.3   1647.8   -790.1   1580.3      173 
## 
## Random effects:
## 
## Conditional model:
##  Groups   Name             Variance Std.Dev. Corr        
##  FISHID   (Intercept)       369.934 19.234               
##           TRIALHypoxia      341.606 18.483   -0.54       
##           TRIALLowSalinity 1032.463 32.132   -0.17 -0.06 
##  Residual                     8.349  2.889               
## Number of obs: 180, groups:  FISHID, 60
## 
## Dispersion estimate for gaussian family (sigma^2): 8.35 
## 
## Conditional model:
##                                                  Estimate Std. Error z value
## (Intercept)                                        52.290      2.511  20.825
## TRIALHypoxia                                      -56.118      2.444 -22.964
## TRIALLowSalinity                                  -42.055      4.182 -10.057
## scale(SMR_contr, scale = FALSE)                   -19.840      4.022  -4.933
## TRIALHypoxia:scale(SMR_contr, scale = FALSE)       11.770      3.914   3.007
## TRIALLowSalinity:scale(SMR_contr, scale = FALSE)   10.884      6.698   1.625
##                                                  Pr(>|z|)    
## (Intercept)                                       < 2e-16 ***
## TRIALHypoxia                                      < 2e-16 ***
## TRIALLowSalinity                                  < 2e-16 ***
## scale(SMR_contr, scale = FALSE)                  8.08e-07 ***
## TRIALHypoxia:scale(SMR_contr, scale = FALSE)      0.00264 ** 
## TRIALLowSalinity:scale(SMR_contr, scale = FALSE)  0.10414    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

cov2cor(vcov(norin.glmmTMB2c)$cond)

##                                                    (Intercept)  TRIALHypoxia
## (Intercept)                                       1.000000e+00 -5.472502e-01
## TRIALHypoxia                                     -5.472502e-01  1.000000e+00
## TRIALLowSalinity                                 -1.759538e-01 -3.991425e-02
## scale(SMR_contr, scale = FALSE)                   3.764461e-15 -1.688551e-15
## TRIALHypoxia:scale(SMR_contr, scale = FALSE)     -9.655529e-16 -4.221495e-18
## TRIALLowSalinity:scale(SMR_contr, scale = FALSE) -8.500783e-16  6.952161e-16
##                                                  TRIALLowSalinity
## (Intercept)                                         -1.759538e-01
## TRIALHypoxia                                        -3.991425e-02
## TRIALLowSalinity                                     1.000000e+00
## scale(SMR_contr, scale = FALSE)                     -1.161186e-15
## TRIALHypoxia:scale(SMR_contr, scale = FALSE)         7.499722e-16
## TRIALLowSalinity:scale(SMR_contr, scale = FALSE)    -2.198129e-15
##                                                  scale(SMR_contr, scale = FALSE)
## (Intercept)                                                         2.642206e-15
## TRIALHypoxia                                                       -1.904170e-15
## TRIALLowSalinity                                                   -1.350437e-15
## scale(SMR_contr, scale = FALSE)                                     1.000000e+00
## TRIALHypoxia:scale(SMR_contr, scale = FALSE)                       -5.472502e-01
## TRIALLowSalinity:scale(SMR_contr, scale = FALSE)                   -1.759538e-01
##                                                  TRIALHypoxia:scale(SMR_contr, scale = FALSE)
## (Intercept)                                                                     -1.892885e-15
## TRIALHypoxia                                                                     4.931878e-16
## TRIALLowSalinity                                                                 1.218256e-15
## scale(SMR_contr, scale = FALSE)                                                 -5.472502e-01
## TRIALHypoxia:scale(SMR_contr, scale = FALSE)                                     1.000000e+00
## TRIALLowSalinity:scale(SMR_contr, scale = FALSE)                                -3.991425e-02
##                                                  TRIALLowSalinity:scale(SMR_contr, scale = FALSE)
## (Intercept)                                                                         -1.281649e-15
## TRIALHypoxia                                                                         9.201455e-16
## TRIALLowSalinity                                                                    -3.060103e-15
## scale(SMR_contr, scale = FALSE)                                                     -1.759538e-01
## TRIALHypoxia:scale(SMR_contr, scale = FALSE)                                        -3.991425e-02
## TRIALLowSalinity:scale(SMR_contr, scale = FALSE)                                     1.000000e+00

norin.glmmTMB2c %>%
  vcov() %>%
  `[[`("cond") %>%
  cov2cor() %>%
  round(3)

##                                                  (Intercept) TRIALHypoxia
## (Intercept)                                            1.000       -0.547
## TRIALHypoxia                                          -0.547        1.000
## TRIALLowSalinity                                      -0.176       -0.040
## scale(SMR_contr, scale = FALSE)                        0.000        0.000
## TRIALHypoxia:scale(SMR_contr, scale = FALSE)           0.000        0.000
## TRIALLowSalinity:scale(SMR_contr, scale = FALSE)       0.000        0.000
##                                                  TRIALLowSalinity
## (Intercept)                                                -0.176
## TRIALHypoxia                                               -0.040
## TRIALLowSalinity                                            1.000
## scale(SMR_contr, scale = FALSE)                             0.000
## TRIALHypoxia:scale(SMR_contr, scale = FALSE)                0.000
## TRIALLowSalinity:scale(SMR_contr, scale = FALSE)            0.000
##                                                  scale(SMR_contr, scale = FALSE)
## (Intercept)                                                                0.000
## TRIALHypoxia                                                               0.000
## TRIALLowSalinity                                                           0.000
## scale(SMR_contr, scale = FALSE)                                            1.000
## TRIALHypoxia:scale(SMR_contr, scale = FALSE)                              -0.547
## TRIALLowSalinity:scale(SMR_contr, scale = FALSE)                          -0.176
##                                                  TRIALHypoxia:scale(SMR_contr, scale = FALSE)
## (Intercept)                                                                             0.000
## TRIALHypoxia                                                                            0.000
## TRIALLowSalinity                                                                        0.000
## scale(SMR_contr, scale = FALSE)                                                        -0.547
## TRIALHypoxia:scale(SMR_contr, scale = FALSE)                                            1.000
## TRIALLowSalinity:scale(SMR_contr, scale = FALSE)                                       -0.040
##                                                  TRIALLowSalinity:scale(SMR_contr, scale = FALSE)
## (Intercept)                                                                                 0.000
## TRIALHypoxia                                                                                0.000
## TRIALLowSalinity                                                                            0.000
## scale(SMR_contr, scale = FALSE)                                                            -0.176
## TRIALHypoxia:scale(SMR_contr, scale = FALSE)                                               -0.040
## TRIALLowSalinity:scale(SMR_contr, scale = FALSE)                                            1.000

Conclusions:

• the average change in SMR for the High Temperature trial at the average control SMR is 52.29.
• on average, the change in SMR during Hypoxia is 56.12 less at the average control SMR than during the High Temperature trial.
• on average, the change in SMR during Low Salinity is 42.06 less at the average control SMR than during the High Temperature trial.
• during the High Temperature trial, for every one unit change in control SMR, the change in SMR declined by 19.84
• the rate of change in SMR during the Hypoxia trial was 11.77 more than that associated with High Temperature.
• the rate of change in SMR during the Low Salinity trial was 10.88 more than that associated with High Temperature (although this was not found to be significant).
• there is approximately twice as much variation within Fish (replicate level) as there is between Fish (2.89 vs
  19.23).

confint

norin.glmmTMB2c %>% confint()

##                                                          2.5 %        97.5 %
## (Intercept)                                       4.736902e+01  5.721163e+01
## TRIALHypoxia                                     -6.090730e+01 -5.132813e+01
## TRIALLowSalinity                                 -5.025118e+01 -3.385951e+01
## scale(SMR_contr, scale = FALSE)                  -2.772227e+01 -1.195774e+01
## TRIALHypoxia:scale(SMR_contr, scale = FALSE)      4.099059e+00  1.944164e+01
## TRIALLowSalinity:scale(SMR_contr, scale = FALSE) -2.242694e+00  2.401121e+01
## Std.Dev.(Intercept)                               1.783349e-05  2.074376e+07
## Std.Dev.TRIALHypoxia                              1.589532e-12  2.149100e+14
## Std.Dev.TRIALLowSalinity                          1.525316e-03  6.768844e+05
## Cor.TRIALHypoxia.(Intercept)                     -9.940613e-01  9.919646e-01
## Cor.TRIALLowSalinity.(Intercept)                 -6.845945e-01  6.843812e-01
## Cor.TRIALLowSalinity.TRIALHypoxia                 6.014078e-01  7.708813e-01
##                                                     Estimate
## (Intercept)                                       52.2903273
## TRIALHypoxia                                     -56.1177117
## TRIALLowSalinity                                 -42.0553426
## scale(SMR_contr, scale = FALSE)                  -19.8400009
## TRIALHypoxia:scale(SMR_contr, scale = FALSE)      11.7703509
## TRIALLowSalinity:scale(SMR_contr, scale = FALSE)  10.8842586
## Std.Dev.(Intercept)                               19.2336575
## Std.Dev.TRIALHypoxia                              18.4825955
## Std.Dev.TRIALLowSalinity                          32.1319577
## Cor.TRIALHypoxia.(Intercept)                      -0.5432691
## Cor.TRIALLowSalinity.(Intercept)                  -0.1658526
## Cor.TRIALLowSalinity.TRIALHypoxia                 -0.0552648

tidy

norin.glmmTMB2c %>% tidy()

## There seems to be a bug, doesn't work with random effects here
norin.glmmTMB2c %>% tidy(effects = "fixed", conf.int = TRUE)

norin.glmmTMB2c %>%
  tidy(effects = "fixed", conf.int = TRUE) %>%
  kable()

effect	component	term	estimate	std.error	statistic	p.value	conf.low	conf.high
fixed	cond	(Intercept)	52.29033	2.510916	20.825196	0.0000000	47.369022	57.21163
fixed	cond	TRIALHypoxia	-56.11771	2.443711	-22.964138	0.0000000	-60.907297	-51.32813
fixed	cond	TRIALLowSalinity	-42.05534	4.181626	-10.057175	0.0000000	-50.251179	-33.85951
fixed	cond	scale(SMR_contr, scale = FALSE)	-19.84000	4.021637	-4.933314	0.0000008	-27.722265	-11.95774
fixed	cond	TRIALHypoxia:scale(SMR_contr, scale = FALSE)	11.77035	3.913997	3.007246	0.0026363	4.099058	19.44164
fixed	cond	TRIALLowSalinity:scale(SMR_contr, scale = FALSE)	10.88426	6.697548	1.625111	0.1041389	-2.242694	24.01121

Conclusions:

• the average change in SMR for the High Temperature trial at the average control SMR is 52.29.
• on average, the change in SMR during Hypoxia is 56.12 less at the average control SMR than during the High Temperature trial.
• on average, the change in SMR during Low Salinity is 42.06 less at the average control SMR than during the High Temperature trial.
• during the High Temperature trial, for every one unit change in control SMR, the change in SMR declined by 19.84
• the rate of change in SMR during the Hypoxia trial was 11.77 more than that associated with High Temperature.
• the rate of change in SMR during the Low Salinity trial was 10.88 more than that associated with High Temperature (although this was not found to be significant).
• there is approximately twice as much variation within Fish (replicate level) as there is between Fish (2.89 vs
  19.23).

tab_model

# warning this is only appropriate for html output
norin.glmmTMB2c %>% sjPlot::tab_model(show.se = TRUE, show.aic = TRUE)

	CHANGE
Predictors	Estimates	std. Error	CI	p
(Intercept)	52.29	2.51	47.37 – 57.21	<0.001
TRIAL [Hypoxia]	-56.12	2.44	-60.91 – -51.33	<0.001
TRIAL [LowSalinity]	-42.06	4.18	-50.25 – -33.86	<0.001
SMR contr	-19.84	4.02	-27.72 – -11.96	<0.001
TRIAL [Hypoxia] * SMR contr	11.77	3.91	4.10 – 19.44	0.003
TRIAL [LowSalinity] * SMR contr	10.88	6.70	-2.24 – 24.01	0.104
Random Effects
σ2	8.35
τ00 FISHID	369.93
τ11 FISHID.TRIALHypoxia	341.61
τ11 FISHID.TRIALLowSalinity	1032.46
ρ01	-0.54
	-0.17
ICC	0.99
N FISHID	60
Observations	180
Marginal R2 / Conditional R2	0.501 / 0.993
AIC	1606.275

Predictions / further analyses

The presence of a significant interaction suggests that the nature of the relationship between SMR change and the control SMR differs between the three trials (High Temperature, Hypoxia and Low Salinity). Equivalently, the differences between the three trials is not constant over all levels of control SMR.

To further tease these interactions apart, we could either:

• compare the slopes separately for each of the three trials
• compare the effects of the three trials at multiple values of control SMR

We will now explore each of the above, although typically we would only do one or the other (depending on which variable was of more interest).

lme (nlme)

Compare slopes

norin.lme2c %>%
  emtrends(~TRIAL, var = "SMR_contr") %>%
  pairs()

##  contrast                      estimate   SE  df t.ratio p.value
##  HighTemperature - Hypoxia       -13.51 3.99 116  -3.384  0.0028
##  HighTemperature - LowSalinity   -10.19 6.76 116  -1.508  0.2911
##  Hypoxia - LowSalinity             3.32 7.97 116   0.417  0.9088
## 
## Degrees-of-freedom method: containment 
## P value adjustment: tukey method for comparing a family of 3 estimates

Conclusions:

• the partial slope associated with control SMR during High Temperature was found to be 13.51 units more negative than during Hypoxia (although not significant). Note, we already had this partial slope in the original summary table.
• the partial slope associated with control SMR during High Temperature was found to be 10.19 units more negative than during Low Salinity (although not significant). Note, we already had this partial slope in the original summary table.
• the partial slope associated with control SMR during Hypoxia was found to be 3.32 units less negative than during Low Salinity (although not significant). Note, we already had this partial slope in the original summary table.

Compare Trials (at different values of control SMR)

We will estimated the pairwise differences between each of the three trials at the minimum, mean and maximum control SMR.

norin.grid <- with(norin, list(SMR_contr = Hmisc::smean.sdl(SMR_contr)))
norin.grid

## $SMR_contr
##     Mean    Lower    Upper 
## 5.178857 3.926671 6.431044

norin.lme2c %>%
  emmeans(~ TRIAL | SMR_contr, at = norin.grid) %>%
  pairs()

## SMR_contr = 5.18:
##  contrast                      estimate    SE  df t.ratio p.value
##  HighTemperature - Hypoxia        50.23  2.49 116  20.150  <.0001
##  HighTemperature - LowSalinity    42.22  4.22 116  10.006  <.0001
##  Hypoxia - LowSalinity            -8.00  4.98 116  -1.608  0.2462
## 
## SMR_contr = 3.93:
##  contrast                      estimate    SE  df t.ratio p.value
##  HighTemperature - Hypoxia        67.14  5.59 116  12.020  <.0001
##  HighTemperature - LowSalinity    54.98  9.46 116   5.814  <.0001
##  Hypoxia - LowSalinity           -12.16 11.15 116  -1.091  0.5217
## 
## SMR_contr = 6.43:
##  contrast                      estimate    SE  df t.ratio p.value
##  HighTemperature - Hypoxia        33.31  5.59 116   5.963  <.0001
##  HighTemperature - LowSalinity    29.46  9.46 116   3.116  0.0065
##  Hypoxia - LowSalinity            -3.84 11.15 116  -0.345  0.9366
## 
## Degrees-of-freedom method: containment 
## P value adjustment: tukey method for comparing a family of 3 estimates

Conclusions:

• as control SMR increases, the magnitude of differences between the trials becomes reduced.
• the High Temperature trial resulted in the greatest positive change in SMR.
• the change in SMR was not found to be different between the Hypoxia and Low Salinity trials.

\(R^2\)

norin.lme2c %>% r.squaredGLMM()

##            R2m      R2c
## [1,] 0.4942092 0.935456

## Nakagawa's R2
norin.lme2c %>% performance::r2_nakagawa()

## # R2 for Mixed Models
## 
##   Conditional R2: 0.935
##      Marginal R2: 0.494

lmer (lmer)

Compare slopes

norin.lmer2c %>%
  emtrends(~TRIAL, var = "SMR_contr") %>%
  pairs()

##  contrast                      estimate   SE df t.ratio p.value
##  HighTemperature - Hypoxia       -13.51 3.99 58  -3.384  0.0036
##  HighTemperature - LowSalinity   -10.19 6.76 58  -1.508  0.2949
##  Hypoxia - LowSalinity             3.32 7.97 58   0.417  0.9088
## 
## Degrees-of-freedom method: kenward-roger 
## P value adjustment: tukey method for comparing a family of 3 estimates

norin.emt <- norin.lmer2c %>%
  emtrends(~TRIAL, var = "SMR_contr") %>%
  pairs() %>%
  as.data.frame()

Conclusions:

• the partial slope associated with control SMR during High Temperature was found to be 13.51 units more negative than during Hypoxia (although not significant). Note, we already had this partial slope in the original summary table.
• the partial slope associated with control SMR during High Temperature was found to be 10.19 units more negative than during Low Salinity (although not significant). Note, we already had this partial slope in the original summary table.
• the partial slope associated with control SMR during Hypoxia was found to be 3.32 units less negative than during Low Salinity (although not significant). Note, we already had this partial slope in the original summary table.

Compare Trials (at different values of control SMR)

norin.grid <- with(norin, list(SMR_contr = Hmisc::smean.sdl(SMR_contr)))
norin.grid

## $SMR_contr
##     Mean    Lower    Upper 
## 5.178857 3.926671 6.431044

norin.lmer2c %>%
  emmeans(~ TRIAL | SMR_contr, at = norin.grid) %>%
  pairs()

## SMR_contr = 5.18:
##  contrast                      estimate    SE df t.ratio p.value
##  HighTemperature - Hypoxia        50.23  2.49 58  20.150  <.0001
##  HighTemperature - LowSalinity    42.22  4.22 58  10.006  <.0001
##  Hypoxia - LowSalinity            -8.00  4.98 58  -1.608  0.2503
## 
## SMR_contr = 3.93:
##  contrast                      estimate    SE df t.ratio p.value
##  HighTemperature - Hypoxia        67.14  5.59 58  12.020  <.0001
##  HighTemperature - LowSalinity    54.98  9.46 58   5.814  <.0001
##  Hypoxia - LowSalinity           -12.16 11.15 58  -1.091  0.5236
## 
## SMR_contr = 6.43:
##  contrast                      estimate    SE df t.ratio p.value
##  HighTemperature - Hypoxia        33.31  5.59 58   5.963  <.0001
##  HighTemperature - LowSalinity    29.46  9.46 58   3.116  0.0079
##  Hypoxia - LowSalinity            -3.84 11.15 58  -0.345  0.9367
## 
## Degrees-of-freedom method: kenward-roger 
## P value adjustment: tukey method for comparing a family of 3 estimates

Conclusions:

• as control SMR increases, the magnitude of differences between the trials becomes reduced.
• the High Temperature trial resulted in the greatest positive change in SMR.
• the change in SMR was not found to be different between the Hypoxia and Low Salinity trials.

\(R^2\)

norin.lmer2c %>% r.squaredGLMM()

##            R2m      R2c
## [1,] 0.4942089 0.876762

## Nakagawa's R2
norin.lmer2c %>% performance::r2_nakagawa()

## # R2 for Mixed Models
## 
##   Conditional R2: 0.877
##      Marginal R2: 0.494

glmmTMB (glmmTMB)

Compare slopes

norin.glmmTMB2c %>%
  emtrends(~TRIAL, var = "SMR_contr") %>%
  pairs() %>%
  summary(infer = TRUE)

norin.emt <- norin.glmmTMB2c %>%
  emtrends(~TRIAL, var = "SMR_contr") %>%
  pairs() %>%
  as.data.frame()

Conclusions:

• the partial slope associated with control SMR during High Temperature was found to be 11.77 units more negative than during Hypoxia (although not significant). Note, we already had this partial slope in the original summary table.
• the partial slope associated with control SMR during High Temperature was found to be 10.88 units more negative than during Low Salinity (although not significant). Note, we already had this partial slope in the original summary table.
• the partial slope associated with control SMR during Hypoxia was found to be 0.89 units less negative than during Low Salinity (although not significant). Note, we already had this partial slope in the original summary table.

Compare Trials (at different values of control SMR)

norin.grid <- with(norin, list(SMR_contr = Hmisc::smean.sdl(SMR_contr)))
norin.grid

## $SMR_contr
##     Mean    Lower    Upper 
## 5.178857 3.926671 6.431044

norin.glmmTMB2c %>%
  emmeans(~ TRIAL | SMR_contr, at = norin.grid) %>%
  pairs() %>%
  summary(infer = TRUE)

Conclusions:

• as control SMR increases, the magnitude of differences between the trials becomes reduced.
• the High Temperature trial resulted in the greatest positive change in SMR.
• the change in SMR was not found to be different between the Hypoxia and Low Salinity trials.

\(R^2\)

norin.glmmTMB2c %>% r.squaredGLMM()

##           R2m       R2c
## [1,] 0.501083 0.9934839

## Nakagawa's R2
norin.glmmTMB2c %>% performance::r2_nakagawa()

## # R2 for Mixed Models
## 
##   Conditional R2: 0.993
##      Marginal R2: 0.501

Summary figures

lme (nlme)

norin.grid <- with(norin, list(SMR_contr = modelr::seq_range(SMR_contr, n = 100)))
newdata <- norin.lme2c %>%
  emmeans(~ SMR_contr | TRIAL, at = norin.grid) %>%
  as.data.frame()
head(newdata)

ggplot(data = newdata, aes(y = emmean, x = SMR_contr)) +
  geom_ribbon(aes(ymin = lower.CL, ymax = upper.CL, fill = TRIAL), alpha = 0.3) +
  geom_line(aes(, color = TRIAL)) +
  theme_classic() +
  theme(
    legend.position = c(0.99, 0.99),
    legend.justification = c(1, 1)
  )

## The .fixed values are the predicted values without random effects
obs <- norin.lme2c %>%
  augment() %>%
  mutate(PartialObs = .fixed + .resid)

ggplot(data = newdata, aes(y = emmean, x = SMR_contr)) +
  geom_ribbon(aes(ymin = lower.CL, ymax = upper.CL, fill = TRIAL), alpha = 0.3) +
  geom_line(aes(, color = TRIAL)) +
  geom_point(data = obs, aes(y = PartialObs, color = TRIAL)) +
  geom_point(data = norin, aes(y = CHANGE), color = "gray") +
  theme_classic()

g1 <- ggplot(data = newdata, aes(y = emmean, x = SMR_contr)) +
  geom_ribbon(aes(ymin = lower.CL, ymax = upper.CL, fill = TRIAL), alpha = 0.3) +
  geom_line(aes(, color = TRIAL)) +
  geom_point(data = obs, aes(y = PartialObs, color = TRIAL)) +
  theme_classic()
g2 <- ggplot(data = newdata, aes(y = emmean, x = SMR_contr)) +
  geom_ribbon(aes(ymin = lower.CL, ymax = upper.CL, fill = TRIAL), alpha = 0.3) +
  geom_line(aes(, color = TRIAL)) +
  geom_point(data = norin, aes(y = CHANGE, color = TRIAL)) +
  theme_classic()
g1 + g2

lmer (lme4)

norin.grid <- with(norin, list(SMR_contr = modelr::seq_range(SMR_contr, n = 100)))
newdata <- norin.lmer2c %>%
  emmeans(~ SMR_contr | TRIAL, at = norin.grid) %>%
  as.data.frame()
head(newdata)

ggplot(data = newdata, aes(y = emmean, x = SMR_contr)) +
  geom_ribbon(aes(ymin = lower.CL, ymax = upper.CL, fill = TRIAL), alpha = 0.3) +
  geom_line(aes(, color = TRIAL)) +
  theme_classic() +
  theme(
    legend.position = c(0.99, 0.99),
    legend.justification = c(1, 1)
  )

## The .fixed values are the predicted values without random effects
obs <- norin.lmer2c %>%
  augment() %>% ## unfortunately, augment.lmer does not back-transform the scale()
  dplyr::rename(SMR_contr = contains("SMR_contr")) %>%
  mutate(
    PartialObs = .fixed + .resid,
    unscale(SMR_contr),
    SMR_contr = V1
  ) # because ggfortify::unscale() returns a data.frame

ggplot(data = newdata, aes(y = emmean, x = SMR_contr)) +
  geom_ribbon(aes(ymin = lower.CL, ymax = upper.CL, fill = TRIAL), alpha = 0.3) +
  geom_line(aes(, color = TRIAL)) +
  geom_point(data = obs, aes(y = PartialObs, color = TRIAL)) +
  geom_point(data = norin, aes(y = CHANGE), color = "gray") +
  theme_classic()

g1 <- ggplot(data = newdata, aes(y = emmean, x = SMR_contr)) +
  geom_ribbon(aes(ymin = lower.CL, ymax = upper.CL, fill = TRIAL), alpha = 0.3) +
  geom_line(aes(, color = TRIAL)) +
  geom_point(data = obs, aes(y = PartialObs, color = TRIAL)) +
  theme_classic()
g2 <- ggplot(data = newdata, aes(y = emmean, x = SMR_contr)) +
  geom_ribbon(aes(ymin = lower.CL, ymax = upper.CL, fill = TRIAL), alpha = 0.3) +
  geom_line(aes(, color = TRIAL)) +
  geom_point(data = norin, aes(y = CHANGE, color = TRIAL)) +
  theme_classic()
g1 + g2

glmmTMB (glmmTMB)

norin.grid <- with(norin, list(SMR_contr = modelr::seq_range(SMR_contr, n = 100)))
newdata <- norin.glmmTMB2c %>%
  emmeans(~ SMR_contr | TRIAL, at = norin.grid) %>%
  as.data.frame()
head(newdata)

ggplot(data = newdata, aes(y = emmean, x = SMR_contr)) +
  geom_ribbon(aes(ymin = lower.CL, ymax = upper.CL, fill = TRIAL), alpha = 0.3) +
  geom_line(aes(, color = TRIAL)) +
  theme_classic() +
  theme(
    legend.position = c(0.99, 0.99),
    legend.justification = c(1, 1)
  )

## The .fixed values are the predicted values without random effects
obs <- norin %>%
  mutate(
    .fixed = predict(norin.glmmTMB2c, re.form = NA),
    .resid = residuals(norin.glmmTMB2c),
    PartialObs = .fixed + .resid
  )

ggplot(data = newdata, aes(y = emmean, x = SMR_contr)) +
  geom_ribbon(aes(ymin = lower.CL, ymax = upper.CL, fill = TRIAL), alpha = 0.3) +
  geom_line(aes(, color = TRIAL)) +
  geom_point(data = obs, aes(y = PartialObs, color = TRIAL)) +
  geom_point(data = norin, aes(y = CHANGE), color = "gray") +
  theme_classic()

g1 <- ggplot(data = newdata, aes(y = emmean, x = SMR_contr)) +
  geom_ribbon(aes(ymin = lower.CL, ymax = upper.CL, fill = TRIAL), alpha = 0.3) +
  geom_line(aes(, color = TRIAL)) +
  geom_point(data = obs, aes(y = PartialObs, color = TRIAL)) +
  theme_classic()
g2 <- ggplot(data = newdata, aes(y = emmean, x = SMR_contr)) +
  geom_ribbon(aes(ymin = lower.CL, ymax = upper.CL, fill = TRIAL), alpha = 0.3) +
  geom_line(aes(, color = TRIAL)) +
  geom_point(data = norin, aes(y = CHANGE, fill = TRIAL), color = "black", shape = 21) +
  theme_classic()
g1 + g2

References

Norin, Tommy, Hans Malte, and Timothy D. Clark. 2015. “Differential Plasticity of Metabolic Rate Phenotypes in a Tropical Fish Facing Environmental Change.” Functional Ecology 30 (3): 369–78. https://doi.org/10.1111/1365-2435.12503.